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Abstract: This paper reports a modeling methodology to predict the effects of operating 
conditions on the thermal behavior of a lithium-ion battery (LIB) module. The potential 
and current density distributions on the electrodes of an LIB cell are predicted as a function 
of discharge time based on the principle of charge conservation. By using the modeling 
results of the potential and current density distributions of the LIB cell, the non-uniform 
distribution of the heat generation rate in a single LIB cell within the module is calculated. 
Based on the heat generation rate in the single LIB cell determined as a function of the 
position on the electrode and time, a three-dimensional thermal modeling of an LIB 
module is performed to calculate the three-dimensional velocity, pressure, and temperature 
distributions within the LIB module as a function of time at various operating conditions. 
Thermal modeling of an LIB module is validated by the comparison between the 
experimental measurements and the modeling results. The effect of the cooling condition 
of the LIB module on the temperature rise of the LIB cells within the module and the 
uniformity of the distribution of the cell temperatures are analyzed quantitatively based on 
the modeling results. 
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1. Introduction 

The lithium-ion battery (LIB) is a preferred power source for hybrid electric vehicle (HEV) 
applications, because it has a relatively higher energy density, longer cycle life, and lower self-discharge 
rate as compared to other cell chemistries. As battery temperature greatly affects the performance, 
safety and life of the LIB for HEV applications, automakers and battery suppliers are paying increased 
attention to the thermal management for an LIB module [1,2]. In a battery module for HEV 
applications, depending on the driving conditions and the type of heating and cooling, an uneven 
temperature distribution can be created in the module. An uneven temperature distribution in an LIB 
module could lead to an electrical imbalance in the battery module and thus lead to lower performance and 
shorter battery life. It is, therefore, important to calculate accurately the uneven temperature distribution in 
a battery module based on modeling to achieve the optimum performance of the battery module [3]. 

There have been many previous works on the modeling of LIBs and the reviews of LIB models are 
given in [4,5]. Doyle et al. [6] developed a porous electrode model of the lithium polymer battery based 
on the concentrated solution theory. Others [7-10] used or modified the model of Doyle et al. [6]. 
Kwon et al. [11] presented a different approach from the rigorous porous electrode model [6-14] to 
predict the discharge behavior of an LIB cell. They developed a two-dimensional model to calculate 
the potential and current density distribution on the electrodes of an LIB cell during constant current 
discharge by solving the equations derived from the principle of charge conservation. Although they 
did not include the potential distribution in the electrolyte phase and the lithium ion transport in their 
model, the model maintains the validity even at high discharge rates and it saves significant 
computational time as compared to the rigorous porous electrode model. Kim et al. [12-14] carried out 
two-dimensional thermal modeling to predict the thermal behaviors of an LIB cell during discharge 
and charge on the basis of the potential and current density distributions obtained by following the 
same procedure of Kwon et al. [11]. They resolved two-dimensionally the heat generation rate as a 
function of position on the electrode of battery cell and state of charge. They found that the 
distributions of the heat generation rate become more non-uniform as the discharge/charge currents 
increase under constant-current discharge/charge mode. They reported good agreement between the 
modeling results and experimental IR (Infrared) measurements. Kim et al. [15] and Yi et al. [16] 
extended their thermal model [12,13] to accommodate the dependence of the discharge behavior on the 
environmental temperature. 

Even though Kim et al. [12-14] reported the non-uniform distribution of the heat generation rate in an 
LIB cell during charge and discharge, the heat generation rate in the battery cell was assumed to be 
uniformly distributed in the many papers related to the thermal modeling for an LIB module or pack 
published until recently [17-22]. In this work, a three-dimensional modeling is carried out to investigate 
the effects of operating conditions on the thermal behavior of an LIB module. The non-uniform 
distribution of the heat generation rate in an LIB cell within the module composed of eight LIB cells is 
calculated by following the modeling approach for a single LIB cell adopted by Kim et al. [12-15]. 
Thermal modeling of an LIB module is validated by the comparison between the experimental 
measurements and modeling results. 
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2. Mathematical Model 


The battery module consisting of eight prismatic pouch-type LIB cells fabricated by Top Battery Co. 
(Ichon, Korea) is modeled in this work. The LIB cell has a nominal capacity of 30 Ah and is composed 
of lithium nickel manganese cobalt oxide (NMC) positive electrodes, graphite negative electrodes, 
and porous separators impregnated with plasticized electrolyte. In the battery module, two cells are 
connected in parallel, composing a single strand, and then four strands are connected in series 
(2P 4S configuration). 

The modeling procedure of a single LIB cell adopted in this work is similar to that of the 
two-dimensional model of Kwon et al. [11]. From the continuity of the current on the electrodes 
during discharge, the following Poisson equations for the potentials in the positive and negative 
electrodes can be derived: 

y2y P = ~ r P J inQ P (1) 

V X=+ r nJ inQ n (2) 


where V P and Vn are the potentials (V) of the positive and negative electrodes, respectively, r P and r n 
are the resistances (Q) of the positive and negative electrodes, respectively, and J is the current density 
(current per unit area (A-m 2 )) transferred through the separator from the negative electrode to the 
positive electrode. Q P and Qn denote the computational domains of the positive and negative 
electrodes, respectively. The resistances, r p and r n , are calculated as described in references [11-13]. 
The relevant boundary conditions for V P are: 



dn 


on T 


V n = 0 on T 

n P2 


(3) 

(4) 


where dl dn denotes the gradient in the direction of the outward normal to the boundary. The first 
boundary condition (3) implies that there is no current flow through the boundary (T P i) of the electrode 
other than the tab. The second boundary condition (4) means that the linear current density through the 
tab (T P 2 ) of length L (cm) has a value of I/L. I is the current (A) during discharge. The boundary 
conditions for Vn are: 


av -=o 

dn 

on T 

1 SV,_i 

on T, 


r p dn L 


(5) 

( 6 ) 


The first boundary condition (5) implies the same as in the case of V P . The second boundary 
condition (6) means that the potential at the tab of the negative electrode has a fixed value of zero as 
the reference potential. The solutions to Equations (1) and (2) subject to the associated boundary 
conditions of Equations (3)-(6) are obtained by using a finite element method. The current density (J) 
of Equations (1) and (2) depends on the polarization characteristics of the electrodes. In this work, 
the following polarization expression used by Tiedemann and Newman [23,24] is adopted: 
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J = Y(V p -V n -U ) (7) 

where Y and U are the fitting parameters. As suggested by Gu [25], U and Y are expressed as functions 
of the depth of discharge (DOD) as follows: 

U= a 0 +a, (DOD)+a 2 (DOD) 2 +a 3 (DOD) 3 ^ (8) 

Y = a 4 +a 5 (DOD)+a 6 (DOD) 2 +a 7 (DOD) 3 +a 8 (DOD) 4 +a 9 (DOD) 5 ( 9 ) 

where ao-a 9 are the constants that provide the best fit of the modeling results to the experimental data 
at an environmental temperature. The values of the parameters ao-a 9 used in Equations (8) and (9) are 
given in Table 1. By solving the equations listed above, the distribution of the current density, J, on the 
electrodes can be obtained as a function of the position on the electrode and time. Therefore, the DOD 
varies along with the position on the electrode and time elapsed during discharge. The distribution of 
DOD on the electrode during discharge can be calculated from the distribution of J as follows: 

[ Jdt 

DOD=DOD n + - (10) 

e T 


where DODo is the initial value of DOD when the discharge process begins, t is the discharge time (s), 
and Qt is the theoretical capacity per unit area (Ahm -2 ) of the electrodes. 

After the distributions of the potential and current density on the electrodes of the battery cell are 
obtained, the heat generation rate within the cell, q, can be calculated as a function of the position on 
the electrode and time as follows: 


q = a J 


V -V -T 

' oc ’cell 1 


dV 

OC 

dT 


VV 


+ CL n 


r P 


+ d r 


Kf 


(in 


where a is the specific area of the battery (m _1 ), J is the current density (A-m“ 2 ) calculated by 
Equation (7), Voc is the open-circuit potential of the cell (V), Vceii is the cell voltage (V) obtained 
through the solutions of Equations (1) and (2), a P and a n are the specific area of the positive and 
negative electrodes (in '), respectively, and |W p | and |W n | are the magnitudes of the gradient vectors 

of VV p and W n obtained by Equations (1) and (2) (Am 1 ), respectively. The detailed definition of 

each term on the right-hand side of Equation (11) is described in references [12,13]. 

Once the heat generation rate within the single LIB cell is determined as a function of the position 
on the electrode and time, the information on the heat generation rate of the cell is transferred to the 
computational fluid dynamics software to perform the three-dimensional thermal modeling of an LIB 
module. The three-dimensional velocity, pressure, and temperature distributions within the LIB 
module were calculated as a function of time at various operating conditions by using the 
computational fluid dynamics (CFD) software ANSYS Fluent 14. A two equations high-Reynolds 
number turbulence model (k- s) was utilized. Mass flow inlet and pressure outlet boundary conditions 
are used for the CFD calculations. The inlet boundary conditions for the flow field are set to 0.005781, 
0.01156, 0.01734, and 0.02312 kg-s _1 for the volumetric flow rates of 10, 20, 30, and 40 ft 3 -min _1 . 



Energies 2014, 7 


7590 


respectively. The outlet boundary condition for the flow field is fixed at 1 atm. The inlet and outlet 
boundary conditions for the temperature field are the fixed air temperature of 25 °C and the free 
boundary condition, respectively. The natural-convection boundary condition is used for the outside 
wall of the module in contact with the environmental air and the value of the heat transfer coefficient 
used for calculation is 8 W-m 2 -°C -1 . The finite volume mesh used for calculation has 3,286,950 cells 
and 693,921 nodes. 

Table 1. Values of the fitting parameters ao-a 9 used for calculations. 


Parameter 

Value 

ao(V) 

4.167 

ai (V) 

-1.283 

a 2 (V) 

1.579 

a 3 (V) 

-1.115 

a 4 (Am -2 ) 

47.042 

a 5 (Am -2 ) 

-179.25 

a6 (Am -2 ) 

839.61 

a 7 (Am -2 ) 

-1759.67 

as (Am -2 ) 

1713.13 

a? (Am -2 ) 

-635.27 


3. Results and Discussion 

In order to test the validity of the modeling, the calculated discharge curves based on modeling for 
the 30 Ah battery fabricated by Top Battery Co. are compared with the experimental data on Figure 1. 
The experiments were performed at room temperature. At various discharge rates from 1 to 3 C, 
the experimental discharge curves are in good agreement with the modeling results. After obtaining the 
distributions of the potential and current density on the electrodes during discharge, the temperature 
distributions of the battery can be calculated as a function of time for various discharge rates by using 
the heat generation rate calculated with Equation (11). The temperature distributions based on the 
experimental IR image and the modeling at the discharge times of 450, 900, and 1800 s with 2 C rate 
and at the discharge times of 300, 600, and 1200 s with 3 C rate are shown in Figure 2a-c, and 
Figure 3a-c, respectively. The overall appearances of the temperature distributions from the experiment 
and modeling are in good agreement. It is observed that the temperature near the current collecting tab 
of the positive electrode is higher than that of the negative electrode. This phenomenon is due to the 
fact that the electrical conductivity of the active material of the positive electrode is much lower than 
that of the negative electrode, although both of the current flows near the tabs of the positive and 
negative electrodes are similarly high as discussed in references [12,13]. The maximum, average, and 
minimum temperatures from the experimental measurement for the LIB cell for the discharge rates of 
2 and 3 C are compared with those predicted by the modeling in Figure 4a,b, respectively. The maximum, 
average, and minimum temperatures from the experiment and modeling are in good agreement for the 
whole range of DOD at both of the discharge rates. 
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Figure 1. Comparison between experimental and modeling discharge curves at discharge rates 
of 1, 2, and 3 C. Open symbols are experimental data and solid lines are modeling results. 



Figure 2. Temperature distributions based on the experimental IR (Infrared) image 
(upper) and the modeling (lower) for the LIB cell at the discharge times of (a) 450 s, 
(b) 900 s, and (c) 1800 s with 2 C rate. 
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Figure 3. Temperature distributions based on the experimental IR image (upper) and the 
modeling (lower) for the LIB cell at the discharge times of (a) 300 s, (b) 600 s, and 
(c) 1200 s with 3 C rate. 


Temperature I 
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Figure 4. Comparison between the maximum, average, and minimum temperatures from 
the experiment and modeling at the discharge rates of (a) 2 C and (b) 3 C. 
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The battery module consisting of 8 prismatic pouch-type LIB cells modeled in this work is shown in 
Figure 5a. The schematic diagram of the cross section of the module along the line A-A of Figure 5a is 
given in Figure 5b and the schematic diagram of the cross section along the line B-B of Figure 5a, 
the direction of cooling air flow, and the positions for thermocouple measurement are shown in Figure 5c. 
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Figure 5. (a) Photograph of the LIB module; (b) Schematic diagram of the cross section of 
the module along the line A-A of (a); (c) Schematic diagram of the cross section along the 
line B-B of (a), the direction of cooling air flow, and the positions for thermocouple 
measurement; and (d) Finite volume mesh used for CFD (computational fluid dynamics) 
calculations and the direction of cooling air flow. 



(a) 


r-J - 

Cell 8 

--'- 1 

1 

r-J - 

Cell 7 

---- 1 

-'-1 

r . 

Cell 6 

-1-'-1 

c_ 

r-l— 

Cell 5 

l- J 

--'-1 

L_ 

r-r— 

Cell 4 

f—■—1 

-1-'-1 

I 

1—1— 

Cell 3 

l-t_l 

L— 

r-J— 

Cell 2 

1 — 3 

1 

Cell 1 

^=3 


(b) 


z 

a 







0 

0 


Battery | Plastic | | Air t Thermocouple 

C^| Inlet |C^ Outlet 

(C) 
























































































Energies 2014, 7 


7594 


Figure 5. Cont. 



The positions for thermocouple measurement are the centers of the rectangular pouch-type cells 
viewed from the top. Based on the information on the non-uniform distribution of the heat generation 
rate within the single LIB cell as a function of the position on the electrode and time, the 
three-dimensional velocity, pressure, and temperature distributions within the LIB module were 
calculated as a function of time at various operating conditions by using the computational fluid 
dynamics software ANSYS Fluent 14. The thermal properties of the LIB module compartments used 
for three-dimensional modeling are given in Table 2. 

Table 2. Thermal properties of the LIB module compartments used for three-dimensional modeling. 


Thermal Property (unit) 

Battery 

Aluminum 

Air 

Plastic 

p(kg-m“ 3 ) 

2765 

2705 

1.225 

1360 

Q (J-kg _1 -°c _1 ) 

1394 

900 

1006 

1500 

k x (W-m _1 -°C _1 ) 

25.7 

231 

0.0242 

0.35 

k y (W-m _1 -°C _1 ) 

25.7 

231 

0.0242 

0.35 

k z (W-m 1 -°C 1 ) 

0.794 

231 

0.0242 

0.35 


The cell temperatures measured with thermocouple are compared with those calculated from the 
modeling in Figure 6a,b at various discharge times with the discharge rates of 2 and 3 C, respectively, 
when the module was cooled in the air by natural convection. Upon the progress of discharge, the cell 
temperatures increase as shown in Figure 6a,b. Two battery cells located at the bottom and top of the 
module (#1 and #8) exhibit lower temperature rises than the other cells. The cell temperatures from 
the modeling agree well with those from the experimental measurement for most of the cells except 
the cell at the bottom (#1) at the discharge rates of 2 and 3 C. The discrepancy between the cell 
temperatures from the experiment and modeling for the cell #1, however, increases as the discharge 
proceeds. This may be due to the fact that the convective boundary condition used for modeling does 
not reflect the convective heat transfer to which the bottom of the module is exposed in the 
experiments. In Figures 7 and 8, the temperature distributions of the LIB cells within the module at a 
few cross sections along the y direction (refer Figure 2a for the y direction) are shown at the discharge 
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times of 1680 s with 2 C rate and 960 s with 2 C rate, respectively, for natural convection cooling. 
As shown in Figure 6a,b, the temperatures of the two cells located at the bottom and top of the module 
are lower than those of the other cells. Along the y direction, the cell temperatures at the cross section 
close to the left wall of the module shown in Figure 5b are lower than those at the cross section close 
to the right wall. That seems to be due to the relatively shorter path for heat transfer from the 
atmospheric air surrounding the module to the cross section close to the left wall as compared to that 
close to the right wall as can be seen in Figure 5b. 


Figure 6. Comparison between the cell temperatures within the module from the 
experiment and modeling at various discharge times with the discharge rates of (a) 2 C and 
(b) 3 C for natural convection cooling. 



(a) (b) 


Figure 7. Temperature distributions of the LIB cells within the module at a few cross 
sections along the y direction at the discharge time of 1680 s with the discharge rate of 2 C 
for natural convection cooling. 
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Figure 8. Temperature distributions of the LIB cells within the module at a few cross 
sections along the y direction at the discharge time of 960 s with the discharge rate of 3 C 
for natural convection cooling 
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In Figure 9a,b, the cell temperatures measured with thermocouple are compared with those 
calculated from the modeling at different discharge times with the discharge rates of 2 and 3 C, 
respectively, for various cooling conditions. As compared to the temperature rise for the case of natural 
convection cooling, it was shown that the temperature rise can be lowered more or less 10 °C near the 
end of discharge process with the discharge rates of 2 and 3 C by increasing the volumetric flow rate of 
cooling air up to 40 ft 3 -min -1 , when the module was cooled by forced convection. The temperature 
rises of the two battery cells located at the bottom and top of the module (#1 and #8) are lower than 
those of the other cells even for the forced convection cooling. That effect, however, diminishes as the 
volumetric flow rate of cooling air increases, because the effect of forced convection becomes more 
dominant at a high volumetric flow rate of cooling air than that of natural convection. The cell 
temperatures from the modeling agree well with those from the experimental measurement at the 
discharge rates of 2 and 3 C. The discrepancy between the cell temperatures from the experiment and 
modeling for the cell #1 that appeared for the natural convection cooling disappears for the forced 
convection cooling, because the relevance of the convective boundary condition used for modeling at 
the bottom of the module becomes less important when the effect of forced convection is more 
dominant at a high volumetric flow rate of cooling air than that of natural convection. 

In Figures 10 and 11, the temperature distributions of the LIB cells within the module at a few cross 
sections along the y direction are shown at the discharge times of 1680 s with 2 C rate and 960 s with 
2 C rate, respectively, for various cooling conditions of forced convection. As shown in Figure 9a,b, 
the temperatures of the two cells located at the bottom and top of the module are lower than those of 
the other cells, but that effect becomes less significant as the volumetric flow rate of cooling air 
increases. Along the y direction, the cell temperatures at the cross section close to the left wall of the 
module are lower than those at the cross section close to the right wall. That effect also becomes less 
significant as the volumetric flow rate of cooling air increases, because the effect of forced convection 
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is more dominant at a high volumetric flow rate of cooling air than that of natural convection through 
the side walls of the module. From the results and discussions on Figures 9-11, we can deduce that the 
uniformity of the distribution of the cell temperatures within the module increase as the volumetric 
flow rate of cooling air for forced convection increases. To quantify this effect, the distributions of the 
volume fraction of the LIB cells within the module as a function of temperature at the discharge times 
of 1680 s with 2 C rate and 960 s with 2 C rate in Figures 12 and 13, respectively, for various cooling 
conditions. As shown in Figures 12 and 13, the mean temperatures of the cells within the module 
decrease and the distributions of the cell temperatures become less dispersed by increasing the volumetric 
flow rate of cooling air. From the results presented in Figures 12 and 13, we can analyze quantitatively the 
effect of cooling condition on the temperature rise of the LIB cells within the module and the uniformity of 
the distribution of the cell temperatures. This analysis is important to the optimal design of the thermal 
management system for HEV applications, because the LIB temperature and the distribution uniformity of 
the cell temperature greatly affect the performance, safety and life of the LIB module. 


Figure 9. Comparison between the cell temperatures within the module from the 
experiment and modeling at different discharge times with the discharge rates of (a) 2 C 
and (b) 3 C for various cooling conditions. 



(a) (b) 

Figure 10. Temperature distributions of the LIB cells within the module at the discharge 
time of 1680 s with the discharge rate of 2 C for the forced convection cooling with the 
volumetric flow rates of (a) 10 ft 3 -min -1 , (b) 20 ft 3 -min -1 , (c) 30 ft 3 -min ', and (d) 40 ft 3 -min -1 . 
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Figure 11. Temperature distributions of the LIB cells within the module at the discharge time 
of 960 s with the discharge rate of 3 C for the forced convection cooling with the volumetric 
flow rates of (a) 10 fit 3 -min -1 , (b) 20 ft 3 -min -1 , (c) 30 ft 3 -min -1 , and (d) 40 ft 3 -min -1 . 



Figure 12. Distributions of the volume fraction of the LIB cells within the module as a 
function of temperature at the discharge time of 1680 s with the discharge rate of 2 C for 
various cooling conditions. 
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Figure 13. Distributions of the volume fraction of the LIB cells within the module as a 
function of temperature at the discharge time of 960 s with the discharge rate of 3 C for 
various cooling conditions. 



4. Conclusions 

A modeling procedure is developed to study the effects of operating conditions on the thermal behavior 
of an LIB module. The potential and current density distribution on the electrodes of an LIB cell are 
predicted as a function of discharge time based on the principle of charge conservation. By using the 
modeling results of the potential and current density distributions of the LIB cell, the non-uniform 
distribution of the heat generation rate in a single LIB cell within the module is calculated. 
The temperature distributions from the modeling of a single LIB cell are in good agreement with those 
from the experimental IR measurement at the discharge rates of 2 and 3 C. Based on the heat 
generation rate within the single LIB cell determined as a function of the position on the electrode and 
time, a three-dimensional thermal modeling of an LIB module is performed to calculate the 
three-dimensional velocity, pressure, and temperature distributions within the LIB module as a 
function of time at various operating conditions. Thermal modeling of an LIB module is validated by 
the comparison between the experimental measurements with thermocouple and the modeling results 
during discharge with the discharge rates of 2 and 3 C. The thermal modeling results for the case of 
natural convection cooling of the module are compared to those for the case of forced convection 
cooling at various volumetric flow rate of cooling air to investigate quantitatively the effect of the 
cooling condition of the LIB module on the temperature rise of the LIB cells within the module and the 
uniformity of the distribution of the cell temperatures. The modeling methodology presented in this 
work may contribute to the optimal design of the thermal management system for HEV applications, 
because the LIB temperature and the uniformity of the distribution of the cell temperatures within the 
module are critical to the performance, safety and life of the LIB module. 
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